
stratifiedSampleSize = function(smallCutoff, x, probSmall = .1, 
  probLarge = 0.02, fullResult=F, probExtraLarge=NULL) {


  popCutoff = pmax(0, x-smallCutoff)
  # probExtraLarge is the proportion to sample in the biggest cities, 
  #  if probExtraLarge has n elements, it's the proportions to sample
  #  from the n biggest regions.
  theorder = order(x, decreasing=T)

  popCutoff = popCutoff / pmax(1, popCutoff[theorder[length(probExtraLarge)+1] ])
  
                                        
  probDiff = probSmall - probLarge
  
  propSampled =  pmax(0,probSmall - popCutoff * probDiff)
  
  # assigning probExtraLarge to the biggest regions
  theorder2 = theorder[seq(from=1, len=length(probExtraLarge), by=1)]
  propSampled[theorder2]=probExtraLarge

  sampleSize = x*propSampled

  if(any(diff(sampleSize[theorder])>0))
    warning("sample size isnt increasing with community size,\n decrease cutoff or increase probLarge or decrease probSmall")

  if(fullResult) 
    return(data.frame(size=x, prop=propSampled, sample=sampleSize))
  else 
    return(sum(sampleSize))

}


getCutoff = function(x, probSmall = 0.1, probLarge = 0.02, totalSize = 150000,
  probExtraLarge=NULL) {

  objFun = function(qq) abs(
  stratifiedSampleSize(qq, x, probSmall=probSmall,
    probLarge = probLarge,  fullResult=F,probExtraLarge=probExtraLarge) -
     totalSize)



  
    
       
  theorder=order(x, decreasing=T)  
  maxNotBig = x[theorder[length(probExtraLarge)+1]]  

    cutoff = optimize(objFun, c(0, maxNotBig-1))$minimum   


  result = list(
    cutoff=cutoff,
    sample = round(stratifiedSampleSize(cutoff, x, 
      probSmall=probSmall, probLarge=probLarge, 
      probExtraLarge=probExtraLarge, fullResult=T) )
  )  
  newx=c(seq(min(x), 
        maxNotBig, len=100),
        x[rev(theorder[seq(1,len=length(probExtraLarge),by=1)])])


  result$toplot = 
    stratifiedSampleSize(cutoff, newx, 
      probSmall=probSmall, probLarge=probLarge, 
        probExtraLarge=probExtraLarge, fullResult=T)
  
    
   return(result) 
    

}

.Random.seed=as.integer(c(1,2,3))

mydata=read.table(file="C:/Documents and Settings/hjiang/My Documents/MyProject/cohort_pilot/community selection/targetpop.csv", header=T, sep=",")
 mydata=mydata[1:65,]

populations = mydata[,3]
names(populations) = mydata[,1]
populations = populations[names(populations)!='Timmins']

result = getCutoff(populations, probSmall = 0.08, probLarge = 0.022, 
  probExtraLarge=c(0.02, 0.025, 0.025))

par(mfrow=c(2,1))  
plot(result$toplot$size, result$toplot$prop, xlab="community size", ylab="proportion sampled", type="l")
abline(v=result$cutoff)
plot(result$toplot$size, result$toplot$sample, xlab="community size", ylab="number sampled", type="l")
abline(v=result$cutoff)
sum(result$sample$sample)


populations = c(round(rnorm(4,100000,10000)),
  120000,
  round(runif(5,10000,15000)),
  rep(50000,4),500000,2000000, 5500,
  rep(25000,2),45000,65000,6500,200000,150000,90000,80000)
 temp=stratifiedSampleSize(20000, populations, probExtraLarge=c(0.02,0.025), probLarge=0.04,probSmall = 0.08,fullResult=T)



  


